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In the Bona-Masso formulation, Einstein equations are written as a set of flux conservative first 
order hyperbolic equations that resemble fluid dynamics equations. Based on this formulation, we 
construct a lattice Boltzmann model for Numerical Relativity. Our model is validated with well- 
established tests, showing good agreement with analytical solutions. Furthermore, we show that by 
increasing the relaxation time, we gain stability at the cost of losing accuracy, and by decreasing 
the lattice spacings while keeping a constant numerical diffusivity, the accuracy and stability of our 
simulations improves. Finally, in order to show the potential of our approach a linear scaling law 
for parallelisation with respect to number of CPU cores is demonstrated. Our model represents the 
first step in using lattice kinetic theory to solve gravitational problems. 


I. INTRODUCTION 

General Relativity was introduced by Albert Einstein 
in 1915, and became the first geometric theory of grav¬ 
itation, including 10 non-linear second order differential 
equations. Due to its complexity, few analytical solu¬ 
tions have been found till today, and therefore, numer¬ 
ical methods have played an important role in the last 
decades. Solving the Einstein equations numerically is 
often called Numerical Relativity. 

Many different systems such as gauge waves, gravita¬ 
tional waves and Schwarzschild black hole systems are 
solved using Numerical Relativity dHi. With the 2005 
breakthroughs, the binary Schwarzschild black holes and 
merger simulations last long enough to extract gravita¬ 
tional waves and study their forms ii- Most popu¬ 
lar methods for solving the equations are based on fi¬ 
nite difference schemes combined with mathematical or 
numerical correction techniques 00. While there is 
a big effort in more complicated phenomena, there is 
also much work done in solving numerical problems en¬ 
countered during the numerical simulations. Numerical 
Relativity has mainly two problems besides stability and 
accuracy: insufficient memory, and complicated geome¬ 
tries that lead to emerging singularities 0. Such prob¬ 
lems, especially in highly curved spacetimes, are tackled 
by mathematical corrections, such as conformal and iso¬ 
metric mappings 0, or by numerical methods, such as 
exponentially growing lattice spacings or adaptive mesh 
refinements 0. The complicated geometries and singu¬ 
larities also lead to highly curved spacetimes and some 
simulations are expected to crash at a certain time 0. 

There are several mathematical formulations of Ein¬ 
stein equations 0. In particular, the 3-1-1 decomposi¬ 
tion is the most convenient mathematical model for nu¬ 
merical calculations and it sets the basis for different for¬ 
malisms. Arnowit-Deser-Misner (ADM) and BSSNOK 
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(Baumgarte, Shapiro, Shibata, Nakamura, Oohara, Ko- 
jima) formalisms have proven to be the most stable ones 
at the present [1, 0 0 . In this work, we use the Bona- 
Masso formulation of Einstein equations, since it con¬ 
sists in first order hyperbolic conservation equations for 
the geometric variables 0,[i3 , resembling fluid dynamics 
equations, and therefore, making suitable the use of fluid 
dynamics solvers. 

In this paper, we propose for the first time a lattice 
Boltzmann model to solve Einstein equations. Usually, 
the lattice Boltzmann method allows us to solve the 
Navier-Stokes equations (or any conservation law) to an 
accuracy that depends on the Knudsen number (roughly 
defined as the ratio of lattice spacing to system size) 
[HI. This method has been successfully used to study 
many physical systems using a fraction of computational 
time in comparison with other numerical methods [H- 
Here, we extend the wide applicability of the lat¬ 
tice Boltzmann method to gravitation. Furthermore, we 
investigate the performance of the model with the fol¬ 
lowing well-established tests reported in Ref. [ia. Ex¬ 
pansion of a flat universe: We simulate an expanding 
flat universe and compare the results with the expected 
analytical solution. This example uses an ideal fluid 
energy-momentum tensor and shows that the energy- 
momentum tensor is coupled correctly to the rest of the 
model variables. Gauge wave: We test the ability of 
the scheme to propagate a gauge wave. The gauge wave 
is achieved through a coordinate transformation, which 
does not change the physics at hand. Linear wave: We 
validate the model propagating the amplitude and phase 
of a gravitational wave, i.e. a spatial transverse wave 
propagating in ^-direction and time. Numerical diffusiv¬ 
ity: We investigate the role of numerical diffusivity in our 
simulation errors. While it is normally added artificially 
in other methods through an elliptic equation 0, it is 
present in the lattice Boltzmann model naturally. Paral¬ 
lelization: We show that our lattice Boltzmann model is 
well suited for parallel computing. The improvement in 
the computation duration for a gauge wave simulation is 
also shown. 

Compared to common methods in Numerical Relativ- 
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ity, such as finite differences, lattice Boltzmann methods 
have extra properties that may lead to improvements or 
short comings, which still need to be tested. Some of 
them can be listed as follow: The numerical viscosity in 
Numerical Relativity is an extra term introduced in finite 
difference methods (see Ref. 0) to add numerical stabil¬ 
ity, while it is naturally present in the lattice Boltzmann 
method; our model uses a lattice that offers more isotropy 
and accurate results when it comes to solving hyperbolic 
equations; and finally, for simple cases, the lattice Boltz¬ 
mann method sets time step equal to the lattice spacing, 
i.e. space and time scale linearly. Therefore, the present 
work is just the beginning of a new way of solving Ein¬ 
stein equation and opens up the door for further research 
in Numerical Relativity. 

This paper first gives an introduction to Numerical 
Relativity and the theory of the specific formulation of 
the Einstein equations that we have chosen. The lattice 
Boltzmann method is introduced in the next section. Af¬ 
terwards, the validation tests are described and their re¬ 
sults using our numerical method are analysed. In the 
last section, we summarise our work and discuss possible 
outlooks. 


II. THEORY 

Einstein equations are 10 coupled second order nonlin¬ 
ear differential equations which, in their most compact 
form, are given by 

Ra/3 — ^a/3; (1) 

where Tq/S is the energy-momentum tensor, the four 
metric and Rap the Ricci tensor defined as 

Rcp = - dpK, + ( 2 ) 

with the Christoffel symbols, 

Ta/? = + dpQaX - dxQap)- (3) 

We use Einstein summation convention throughout the 
paper and the greek letters indicate summation over four 
dimensions, while the latin letters indicate summation 
over the three spatial components. Einstein equations, 
Eq. dl]), in their expanded forms possess very complex 
shapes and several formalisms are used to ease their so¬ 
lutions. 


A. 3+1 Arnowit Deser Misner formalism 

The ADM formalism of Einstein equations was first 
introduced through Hamiltonian framework of General 
Relativity [1, [l^ . Here, we obtain the ADM formalism 
with a slightly different construction, as described below. 


The 3 + 1 decomposition of Einstein equations lets us find 
a solution for the four-metric in different hyper-spaces, 
and is given as follows: 

= adt^ — + f5^dt){dx^ + 13^ dt) , (4) 

where a and j3i are describing the evolution of hyper¬ 
spaces in the four dimensional space, and the hyper¬ 
spaces (also called slices) are associated with a certain 
three-metric 7 ^. The construction of these geometric 
objects is given by the following scheme: we construct 
three dimensional spaces with the three-metric 7 ^- and 
find the normal vector to these spaces with components 
a and jdi. At every time step we need a variable that 
describes the evolution of 7 ^^-, which is given by the ex¬ 
trinsic curvature , which is the projection of the cur¬ 
vature to the normal vector onto the three spaces. Due 
to the symmetry of 7 ^ (thus Kij) we arrive to 6 differ¬ 
ential equations describing the evolution of the system. 
Note that other 4 degrees of freedom are still left, as 
we are dealing with 10 differential equations. These de¬ 
grees of freedom are covered by the constraint equations, 
which are energy and momentum conservation. Einally 
the ADM equations and constraints are given by the fol¬ 
lowing equations: 

(dt - Li3)-fij = -2aKij, (5) 

{dt 

a(i?g^ - + tr(A)A,, - ( 6 ) 

- tr{K^) + tr(A)2 - 2a^G°° = 0, (7) 

KG-d,{tT{K))-aG° =0, (8) 

where Lp corresponds to the Lie derivative along the vec¬ 
tor /3 and Gap to the energy momentum tensor, and 
is the projection of the four Ricci tensor on the three di¬ 
mensional spaces. These equations correspond to one 
representation of the Einstein equations. As we will see 
later, there can be other formulations. Indeed, all phys¬ 
ical solutions of the equations match, even if the formu¬ 
lations may not have the same mathematical structure. 
The physical solutions are known to be the solutions be¬ 
longing to the “constrained space”. However, the math¬ 
ematical differences may increase or decrease the stabil¬ 
ity of numerical schemes. These two equations can be 
evolved in two different ways. The first one is letting the 
system evolve freely and monitor the constraints (“free 
evolution”), and the second one, is solving the constraint 
equations for each time step for some of the variables 
and make sure that they are fulfilled (“constrained evo¬ 
lution”). It must be pointed out that the 3+1 decomposi¬ 
tion does not lead to an evolution equation for the slicing 
a and the shift /?,. So we have four additional degrees 
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of freedom, which we can exploit to improve stability in 
different systems. Thus, we can formulate 


{dt - Lp)a = -c?Q, 

(9) 

{dt - Lp)l3i = 

(10) 


where Q and Qi are gauge functions. 

1. Bona-Masso formulation 


equations: 

Slicing: 

dta = a{P^Ar — aQ), (13) 

Three metric: 

9*7,, = - s,,), (14) 


While there are many ways to solve Einstein equations 
numerically, the following two methods are the most pop¬ 
ular ones: ADM and BSSNOK (Baumgarte, Shapiro, 
Shibata, Nakamura, Oohara, Kojiina) formalisms Q. 
While ADM tries to solve the equations directly, BSS¬ 
NOK makes use of conformal mappings and increases 
the stability of the evolution as conformal variables are 
evolved. Both formalisms are solved mainly by finite dif¬ 
ference methods. In this work, we will use a third method 
known as Bona-Masso formalism 0,113) [I3| • It consists 
in a set of first order flux conservative hyperbolic equa¬ 
tions equivalent to Einstein equations. We have seen that 
the Einstein equations, within the ADM formalism, are 
first order in time derivatives but second order in space 
derivatives. Thus, the Bona-Masso formalism introduces 
the following space derivatives to obtain first order dif¬ 
ferential equations: 

Ak = 9fcln(a), Bl = 9fe/372, Dkij = dkHjl^. (11) 

There are two evolution systems in this formalism 0|. 
The first one is the “Ricci evolution system”, which con¬ 
sists in the above formulated ADM equations, and the 
second one is the “Einstein evolution system”. In the 
latter one, the four Ricci tensor is replaced by the energy- 
momentum tensor (using Eq. (UJ). Thus in the Ein¬ 
stein system the energy constraint is already solved by 
the evolution equations. The Bona-Masso formulation 
introduces the following variable similar to the energy 
constraint: 


Efc = DJ - D 7 . ( 12 ) 

The quantity Vi allows us to implement the momen¬ 
tum conservation in its evolution equations such that 
we do not need to solve the momentum constraint sep¬ 
arately, i.e. the momentum conservation is solved by 
the equations for Vi during the evolution. Now we 
can take a look at the Einstein equations through the 
Bona-Masso formalism. As mentioned in the introduc¬ 
tion, we are considering zero shift cases, which means 
/3i = 0, Qi = 0 => Piit) = 0. The expressions include 
the shift vector only for general purposes. 


B. Bona-Masso evolution equations 

We are considering the Einstein evolution system. 
The Bona-Masso equations are given by the following 


Bii -\- Bj 


J1- 


a 


(15) 


Extrinsic Curvature: 

dtK,, + dr{-p''K,, + aXT) = aS,j, (16) 

+ 27 - + 

i5}'(A, + 2E,-D,;)--5^, (17) 

+ trKK,j + 

^ (K^rB/ + K,rB: - IUMB)) + + 2D/^DJ+ 

TlTT - - {20^ - ^7(A,. + A..)+ 

(vj - ^ j + A, (v. - + 

^(-D;*r7 + - 2v^a,,+ 

tr(A:7 - {tvKf + 2a2G'™), (18) 

R'ff = A, - + tr(G)), (19) 

Momentum related variable: 

dtVk + 97-/3”Rfc + B\ - Bn = Dfc = 
a[aGl + Ar{Kl - Xv{K)5l) + (D,; - 2D/^)k: 

-A7(D,|-2D^.;()]+2(Rfc”-5^trR)W+2(D7-<5^Dj7R”, 

(20) 


First order derivatives: 

Ak = dk ln(Q;), Dkij = dklijl2. (21) 

where are the space components of the four¬ 

dimensional Ricci tensor. Although in the original for¬ 
mulation of Bona-Masso equations, Ai and Dijk are con¬ 
sidered independent variables with their own respective 
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evolution equations, we have evolved the slicing and the 
metric first and then taken the derivatives. During the 
simulations we observed that, if they are evolved inde¬ 
pendently, then the derivative of the slicing (the metric) 
does not match with (Dijk) which leads to errors. 
Therefore we will calculate them directly through a and 

lij- 


C. Gauge choices 



The 3-1-1 formulation does not define slice and lapse 
evolution uniquely, which gives us a freedom to choose 
the appropriate slicing and lapse depending on the sys¬ 
tem. The simplest slicing would be the geodesic slicing 
(5 = 0, where the time component of the four metric 
stays constant throughout the simulation. However in 
some situations this slicing fails to support the stability 
of the program, as the evolution of the metric is coupled 
to it. For example, for black holes the geodesic slicing 
leads to a negative metric, thus to a numerically unsta¬ 
ble algorithm @ . Other slicings are: the maximal slicing 
tr(iF) = 0, harmonic slicing Q = tr(iF), and “1 -I- log” 
slicing Q = tr(iF)/a. The maximal and 1 -I- /oo slicings 
are also known as singularity avoiding slicings Q because 
they collapse the time evolution close to the singularity. 
The harmonic slicing corresponds to a sliciM where a 
fulfils the wave equation (9^ — d"^)a = 0 [^. In the 
following tests we only use the geodesic and harmonic 
slicings and work with zero shift {Qi = 0). 

The tests have periodic boundary conditions. We will 
avoid tests that include other boundary conditions, e.g. 
“single Schwarzschild black hole”, since they are known 
to be very challenging and are currently under intense 
research (combined with appropriate slicings and lapses). 


III. LATTICE BOLTZMANN MODEL 

The lattice Boltzmann Method can be used to sim¬ 
ulate fluids or solve partial differential equations in the 
form of conservation laws [2lj. We start with the discrete 
Boltzmann equation for the distribution functions, fx, 

f\{x + CxSt,t + St) - fx{x,t) = ^ 

T/dt 

( 22 ) 

where t is the relaxation time, and the index A repre¬ 
sents the discrete velocity {cAjAeAf- In this work we will 
use the Bhatnagar-Gross-Krook (BGK) approximation 
[HI, which is a small amplitude approximation for the 
collision term of the lattice Boltzmann equation (rhs of 
Eq. (12^ 1. In case of an ideal gas, one can take as 
the Maxwell-Boltzmann (MB) distribution, expanded in 
orthogonal polynomials, and recover the Navier-Stokes 
equations [^. The macroscopic fields can be calculated 


FIG. 1. D2Q9 lattice configuration vectors. The thickness of 
the arrows represent the weights. The dark red point in the 
middle denotes the (0,0, 0) vector. 


by using the relations 

P = Ylf^' = ( 23 ) 

A A 

which, in the case of fluid dynamics, p is the fluid den¬ 
sity and pui the momentum density. The lattices are 
described by their dimensionality (D) and amount of dis¬ 
crete velocity vectors (Q). An n-dimensional lattice with 
m velocities would be denoted by DnQm. Here, we will 
use D2Q9 and D3Q19 in two and three dimensions, re¬ 
spectively. The D2Q9 lattice configuration is given by: 

f \ fO il 0 il il\ 

(Co,Ci_2,C3,4,C5,6,C7,8) = I Q q I (24) 

where the weights are defined by wq = wi_ 2 , 3,4 = 

and the speed of sound Cs = -^ (see figure 
[1]). On the other hand, the D3Q19 lattice configuration 
is given by: 



±1 

0 

0 

±1 

±1 

±1 

±I 

0 

0 \ 


0 

0 

±1 

0 

Tl 

±1 

0 

0 

±1 

±1 

(25) 

VO 

0 

0 

±1 

0 

0 

Tl 

±I 

Tl 

±V 


with weights 

Wo 

= 

W'1.2,3, 

, 4 , 5,6 

= ■ 

IS ’ 

M and 


Cs = Both lattices are accurate up to second order, 
which means that one can recover the moments of the 
distribution up to second order [2^. 

The discrete Boltzmann equation with a source term 
distribution Sx{x,t) is given by 

fx{x + cxSt,t + St)-fx{x,t) =- *~^’(26 ) 

-I- St Sx{x, t). 

This equation describes the evolution of the distribution 
functions. We calculate the equilibrium and source dis¬ 
tributions for each component of a, 7 ij, Vi, and Kij, such 
that the correct moments of the equilibrium distribution 
are satisfied, and consequently, the right macroscopic dif¬ 
ferential equations are recovered. Thus, one gets 

“/r = wxa, 


(27) 
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(28) 


K £eq 

4 a = 

“^A = «^A (^1 - a(/3M, - aQ), (31) 

'^4a = «^a (^2 - (/3"D„, - a{K,j - Si,)), (32) 

""S.A = «^A (^1 - (33) 

^S'ijA = wxaSij ^1 - , (34) 

where the symbol * at the left of the distribution and 
source term, */ and *5, denotes the field to which / and 
S are associated. The macroscopic variables of concern 
are calculated by 

* St*S 

P = /a + (35) 

A 

where p* stands for (a, 7 i,, Vi, Kg). By performing 
the Chapman-Enskog expansion , we see that the 

distributions and source terms recover the Bona-Masso 
formulation of Einstein equations to the first order in 
Knudsen number [U, [2^. One last remark is that we 
could fix the second order moment of the equilibrium 
distributions to zero, but it decreases the stability of the 
system for two reasons: first, the equilibrium distribu¬ 
tion can become easier negative, which violates the H- 
theorem leading to an unstable evolution of the system; 
and second, the numerical viscosity increases the stability 
of the system in the same fashion as the numerical dif- 
fusivity, introduced in other methods, does it. We have 
finished with the model description. Now we will study 
some well-established examples in order to validate and 
characterise the method. 




•riL \ rn 

Ca 


(30) 


Efe 1 - 




(29) 


IV. TESTS AND RESULTS 

In order to validate our model, we will perform three 
numerical tests. The first two tests are flat space-times, 
in one we consider an ideal fluid and in the other one just 
vacuum for the energy-momentum tensor. Finally, the 
last test corresponds to a curved space-time in vacuum. 
The errors of the simulations are presented in section El 


A. Expansion of flat universe 


In this case, the space part of the energy-momentum 
tensor, Gtj, of an ideal fluid is given by 

Goo = p{t)c^, Gij = P(t)jij, (36) 

where p(t) and P{t) are the density and pressure of the 
fluid, respectively. By solving the Einstein equations for 
this system one obtains the Friedmann equations [ 2 ^ 
(also knows as FLRW metric), which are 

-I- kc^ 9 ttGp + h.(? 


a 

a 



Ac2 

“ 3 “’ 


(38) 


where a is a quantity that determines the metric. 


7i,- = a(t)'^Stj = t^Sij, (39) 

and 

pm = T. (40) 

For an ideal gas with equation of state P = p (assuming 
= SttG = 1). The simulation ran using a D3Q19 lattice 
with St = Sr = 0.001, and a relaxation time r/St = 3 (all 
values are given in numerical units). The results can 
be observed in Figs. [2] and |3l We have used a lattice 
of 4 X 4 X 4 cells due to the fact that the lattice size 
does not matter due to isotropy and homogeneity of the 
problem, i.e. the expansion is the same at every position. 
In Figs. [2]and[3l we observe excellent agreement between 
the results of our simulation and the theory, showing that 
the universe expands following the Friedmann equations. 


B. Gauge Wave 

In this test, the metric tensor is given by 

ds^ = gapdx°‘dx^ = —Hdt^ + Hdx^ + dy^ + dz^, (41) 

where H = H{x — t) = 1 — A sin d the 

wavelength of the gauge wave, which in our case is set 
to unity, and A = 10“^ its amplitude. The extrinsic 
curvature can be calculated by directly taking the time 
derivative of the three metric and dividing by —2a, 



It can be checked by Eq. ([61) that the extrinsic curvature 
also evolves with the time derivative of this expression. 
The time evolution of a is given by 

dta = —a^Q Q = tr{K). 


(43) 
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FIG. 2. Time evolution of a single component of the metric 
tensor for an expanding universe. The metric starts from 
unity and evolves according to the time dependence given in 
Eq. (1391) . The lattice size and spacing are irrelevant factors 
due to isotropy and flatness. 



FIG. 3. Single component of the second fundamental form, 
Kxx, for an expanding universe. This fundamental form de¬ 
scribes the time derivative of the metric via K{t) = —a{t)/2. 


Here, we have chosen the harmonic slicing because we 
want a = '/H to propagate also as a wave and the har¬ 
monic slicing provides a wave-like evolution to a, while 
keeping the evolution of the quantity H consistent with 
the evolution of the three metric. If the evolution of a 
and 7 do not lead to the same equations for H, then the 
respective time derivatives of the extrinsic curvature do 
not match. The results, using r/St = 4.5, are presented 
in Figs. m [3 and El We see again that they are in good 
agreement within a relative error of 1%. The simulations 
ran using a D2Q9 lattice configuration, where we have re¬ 


defined the spatial coordinate x ^ Xr and keep x for the 
positions in the lattice, being Xr = —4005r-|-(a; — 1/2) 6r, 
and St = Sr = 0.00125. 



FIG. 4. Evolution of the xx component of the metric for the 
propagation of a gauge wave is shown at three different times. 
Here, t denotes the numerical time step. We see that the 
simulation deviates slightly from the analytical values, but 
the shapes are preserved and the maximum (minimum) of 
the waves are on the same track. “Ana” stands for analytical 
solution and “Sim” for simulation. 



FIG. 5. Evolution of the second fundamental form for the 
propagation of a gauge wave is shown at three different times. 
We see a similar behaviour to that of the metric. “Ana” 
stands for analytical solution and “Sim” for simulation. 


C. Linear Wave 

In this test, we set the following metric tensor: 

ds^ = gai 3 dx°‘dx^ = —dt^ -b dx^ + (1 + d)dy‘^ + (1 — b)dz‘^, 

(44) 
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FIG. 6. Evolution of a {tt component of the metric) for the 
propagation of the gauge wave is shown at three different 
times. Here, we observe again a similar behaviour to that of 
the metric. “Ana” stands for analytical solution and “Sim” 
for simulation. 



FIG. 7. Evolution of the xx component of the metric for the 
propagation of a linear wave is shown at three different times 
t. Note that the linear wave evolves similarly to the gauge 
wave. 


where b = with d the system size as 

above. The extrinsic curvature can be calculated directly 
by taking the time derivative of the three metric and di¬ 
viding by —2a, with a = 1: 


ttA 

^VV ~ 


27r(x — t) 


(45) 


K 


ZZ 



^27r(a; ~ ^ 


(46) 


In this case, we take geodesic slicing <5 = 0, which would 
impose into the system the analytical evolution of a. As 
the perturbation is traceless, geodesic slicing and har¬ 
monic slicing are identical gauges. The amplitude is cho¬ 
sen as A = 10“^. For this simulation, we use a D3Q19 
lattice configuration, r/St = 4.5, and the same values 
for X, t, and St, Sr, as before. Note that the results are 
very similar to those of the gauge wave as expected (see 
Fig. El). 


V. ANALYSIS AND SIMULATION ERRORS 

In this section, we analyse the errors of the performed 
simulations of the previous section. The tests will be 
analysed in groups with the corresponding errors. 


A. Errors in the Validation Tests 

The relative error of the metric for the case of an ex¬ 
panding flat universe is shown in Fig. [51 As it can be 
seen, the error decreases monotonically while the simu¬ 
lation stays in good agreement with the theory. 
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FIG. 8. Relative error of the metric for an expanding universe. 
It is observed that throughout the simulation, it is always 
below 0.01%. 


On the other hand, the gauge and linear waves show 
very similar behaviour in our simulations. The test of 
the gauge wave is expected to get unstable by design 
due to emerging singularities. As mentioned in Ref. [l^ . 
such system with topology must have a singularity 
in the future or in the past and its effect must show at 
some time during the simulation. Our simulation ran for 
almost 2 • 10^ time steps till the relative error exceeds 
10%. Compared to the well-developed methods of ADM 
and BSSNOK with the use of Cactus or other numeri¬ 
cal solvers, which can handle much more time steps, our 
model needs improvements. This error introduced by the 
lattice Boltzmann model is due to the fact that singulari¬ 
ties possess large gradients in the geometric variables and 
they lead to negative equilibrium distribution functions, 
and consequently, to numerical instabilities. Thus, we 
can expect that by adding the H-theorem, i.e. introduc¬ 
ing a entropic lattice Boltzmann model , one can 


















improve drastically the model keeping its simplicity. 


B. Relaxation Time r and Numerical Diffusivity rj 

The relaxation time t is one of the characteristic 
parameters of any single relaxation lattice Boltzmann 
scheme. More precisely, it is related to the numerical 
diffusivity through 



For the above tests, r/St varies between 2.5 and 5.25, 
and consequently, r] takes values between 0.00083 and 
0.002083. We see that as the relaxation time increases, 
the stability of the system also increases (see inset of 
Fig. ini), while the errors become larger (see Fig. HD. For 
instance, we observe that while r/St = 2.5 leads to nu¬ 
merical instabilities quicker than r/St = 4.5, it also intro¬ 
duces a less mean relative error into the system. Indeed, 
further tests show that there is an optimal value of r 
for which the error introduced by the model is less while 
the stability is better. In the case of gravitational waves 
T/6t = 4.5 is the optimal relaxation time. 

Note that Eq. (ITTl) states that by increasing the relax¬ 
ation time, the numerical diffusivity also increases for a 
fixed St = Sx, but one can also fix the numerical diffu¬ 
sivity while changing independently the relaxation time. 
Thus, in a second test (see Fig. [TUD we study the be¬ 
haviour of the error when one increases the relaxation 
time and decreases the lattice spacing, such that the nu¬ 
merical diffusivity remains constant. We observe that the 
simulation lasts about 50% longer in order to obtain the 
same accuracy. Therefore, we conclude that the numer¬ 
ical diffusivity tunes the accuracy, while the relaxation 
time the stability of the model. 

As a final error test, the inset of Fig. [TUI shows that the 
error decreases linearly with the lattice spacing (at time 
t = 400(5t), concluding that our model converges to the 
analytical solution when the lattice spacing is decreased. 


VI. PARALLELISATION OF THE CODE 

Our code has been parallelised with OpenMP and 
tested with 1 — 24 cores. The results are given in the 
Fig. [TTJ The linear decrease in total computational time 
with respect to number of cores shows that our model 
is optimal for parallelisation. Note that for 24 cores, 
the total computational time deviates slightly from the 
linear behaviour because we have performed our simula¬ 
tions in nodes with a maximum of 24 cores, which implies 
that we have reached the physical limit, and usually at 
this limit tasks are not performed very efficiently. Other 
more sophisticated implementations as MPI and CUBA, 
for GPUs, will be a subject of future works. 



FIG. 9. Mean relative error of 'yxx for the propagation of 
gravitational waves at the first 500 time steps. For higher 
relaxation time the accuracy gets worse. However, it can be 
seen, in the inset, that the life-time of the simulation increases 
with increasing the relaxation time. 



FIG. 10. Main frame: relative error of the metric for the 
propagation of gravitational waves at different time steps for 
different resolutions. Here, we have kept constant the numer¬ 
ical diffusivity, 77 = 0.0016. The errors strongly depend on 
the lattice spacing, and finer lattices lead to longer simula¬ 
tion time. Inset: we see that with finer lattices the relative 
error introduced at time step 400 St decreases as a power law 
showing convergency. 

VII. CONCLUSIONS 

In this paper, we have developed a lattice Boltzmann 
model for solving Einstein equations, using the Bona- 
Masso formalism. We have validated our model with 
well-established tests, namely the expansion of a flat uni¬ 
verse, and the propagation of gauge and linear waves. 
The expansion of a flat universe was recovered accurately 
and the wave tests showed good agreement with the an¬ 
alytical solutions. The roles of the relaxation time and 
the numerical diffusivity were also studied hnding that 
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FIG. 11. Computational time as function of the number of 
CPU cores. Here we observe that the computational time falls 
rapidly with increasing the number of CPU cores and is fitted 
very well by an inverse linear function. 

they are crucial in determining the stability and accu¬ 
racy of the model. In particular, we have observed that 


the system gains stability (accuracy) by increasing (de¬ 
creasing) the relaxation time, and therefore, an optimal 
value, that compromises both, can be found. More pre¬ 
cisely, the numerical diffusivity tunes the accuracy while 
the relaxation time the stability of the model. 

In addition to the validation tests, the inverse linear 
dependence of computational time with respect to the 
number of CPU cores was demonstrated, which is a ma¬ 
jor strength of lattice Boltzmann methods. It must be 
clearly underlined that with further work on this model, 
e.g. entropic extensions, the lattice Boltzmann method 
might offer new numerical advantages to Numerical Rel¬ 
ativity. 
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